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Abstract. A mean-field theory is developed for the scale-invariant length dis- 
tributions observed during the coarsening of one-dimensional faceted surfaces. 
This theory closely follows the LSW theory of Ostwald ripening in two-phase 
systems [T] [2] [3] , but the mechanism of coarsening in faceted surfaces requires 
the addition of convolution terms recalling the work of Smoluchowski [4] and 
Schumann [5] on coalescence. The model is solved by the exponential dis- 
tribution, but agreement with experiment is limited by the assumption that 
neighboring facet lengths are uncorrelated. However, the method concisely 
describes the essential processes operating in the scaling state, illuminates a 
clear path for future refinement, and offers a framework for the investigation 
of faceted surfaces evolving under arbitrary dynamics. 



1. Introduction 

In many examples of faceted surface evolution, a facet velocity law giving the 
normal velocity of each facet can be observed, assumed, or derived. Examples of 
such dynamic laws describe growth of polycrystalline diamond films from the vapor 
O [7] , evolution of faceted boundaries between two elastic solids [8] , the evapora- 
tion/condensation mechanism of thermal annealing [5], and various solidification 
systems [TUl EJ HH EH Q3| ■ Such velocity laws are typically configurational, de- 
pending on surface properties of the facet such as area, perimeter, orientation, 
or position, and reduce the computational complexity of evolving a continuous 
surface to the level of a finite-dimensional system of ordinary differential equa- 
tions. This theoretical simplification enables and invites large numerical simula- 
tions for the study of statistical behavior. This has been done frequently for one- 
dimensional surfaces [TDl[niII3II3IIlIIllIlSliniIIHlIISl, while less frequently for 
two-dimensional surfaces due to the necessity of handling complicated topological 
events [20l ETJ [22l [23l [24] . Such inquiries reveal that many of the systems listed 
above exhibit coarsening - the continual vanishing of small facets and the increase 
in the average length of those that remain. Notably, these systems also display dy- 
namic scaling, in which common geometric surface properties approach a constant 
statistical state, which is preserved even as the length scale increases. 

The dynamical scaling behavior of coarsening faceted surfaces recalls the pro- 
cess of Ostwald ripening [25] , in which small solid-phase grains in a liquid matrix 
dissolve, while larger grains accrete the resulting solute and grow. As this pro- 
ceeds, the distribution of relative particle sizes tends to a constant stat^E In fact, 
general 2D faceted surfaces fall conceptually into the same class of phase-ordering 
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systems, except that the a vector order parameter reflecting surface normal re- 
places the scalar order parameter reflecting phase |31j : other systems exhibiting 
similar behavior include coarsening cellular networks describing soap froths and 
polycrystallinc films [32j [33] [34] [35l [36] , and films growing via spiral defect [37] . 
Each of these generalizations is characterized by a network of evolving boundaries 
which separate domains of possibly differing composition, and exhibit coarsening 
and convergence toward scale-invariant steady states. 

Since dynamic scaling pushes complex systems into a state which can be effec- 
tively characterized by just a few statistics, it is natural to seek simplified models 
which replicate this behavior. The canonical example of this approach is the cele- 
brated theory of Lifshitz, Slyozov, and Wagner describing Ostwald ripening [HOE]. 
Generically, such an approach selects a distribution of some quantity, and includes 
just enough of the total system behavior to specify the effective behavior of that 
quantity - for example, the original LSW theory first identifies the average behav- 
ior of particles as a function of size, and uses that result to identify a continuity 
equation describing distribution evolution. Ideas of this kind have been applied 
to several of the higher-order cellular systems introduced above - for soap froths 
[38] [39] , polycrystals [] , and spiral-growth films [37] . To the extent that such ap- 
proaches mirror experimental data, they can yield valuable physical insight which 
cannot be gained by considering single particles, nor even by direct numerical sim- 
ulation of larger ensembles. However, to date no similar attempt has been made 
for evolving faceted interfaces. 

Given the wide variety of examples of purely faceted motion, the membership of 
this problem class in the wider class of phase ordering systems, and the past success 
in applying mean-field analyses to these systems, it is somewhat surprising that no 
such attempt has been made to describe the mean-field evolution of faceted surfaces. 
In this chapter, therefore, we take a first step in that direction by introducing a 
framework for describing the distribution of facet lengths in ID faceted surface 
evolution. Our approach closely resembles the LSW theory of Ostwald ripening, 
in that a facet-velocity law allows the effective behavior of facets by length, and 
thus the specification of a continuity equation governing the evolution of the length 
distribution. However, our model differs in that facets do not vanish in isolation 
as do grains in Ostwald ripening - instead each vanishing facet causes its two 
immediate neighbors to join together. This process of merging is not treated in 
LSW theory, and requires the introduction of a convolution integral reminiscent of 
equations due to Smoluchowski U and Schumann [S] describing coagulation. We 
apply our method to one particular facet dynamics, associated with the directional 
solidification of faceting binary alloys. However, the method is general and can be 
applied to any dynamics where effective facet behavior is accessible. 



2. Example Dynamics and Problem Formulation 

During the directional solidification of a strongly anisotropic binary alloy, small- 
wavelength faceted surfaces develop. If the alloy is solidified above a critical velocity, 
a layer of supercooled liquid is created at the interface, which drives a coarsening 
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instability governed by the facet dynamics [H] 
(1) 



dhi 

It = {kt) 



Figure [T^ displays representative surface evolution during coarsening. There, the 
locations of corners are plotted over time. We see that this system exhibits bi- 
nary coarsening, whereby a single facet shrinks to zero length, causing two corners 
meet and annihilate. As coarsening proceeds, the average facet length increases 
(Figure Ht), and a scale- invariant distribution of relative facet lengths is reached 
(Figure Eh). 




Figure 1. Survey of coarsening behavior, (a) A representative 
example of the kink/anti-kink evolution (red/blue), (b) Facet 
lengthscale growth with time, (c) Scale-invariant distribution of 
relative facet lengths (the tail is gaussian). 

To describe the evolution of a scale-invariant length distribution such as that 
shown in Figure [TJ;, we will think in terms of a mass-transport problem. We will 
derive equations governing the evolution of the mass distribution P{x), where mass 
"particles" at position x correspond to facets with length L = x. Here P(x) is 
not necessarily a probability distribution, and so its mass is not assumed to be 
unity. Thus, it will be useful to also consider the unit-mass scaling of this function, 
P{x) = P{x) J °° P(x)dx, which is a true probability distribution. The goal, 
of course, is to model P, but we will begin with a simple equation governing the 
evolution of P, and add terms until that equation governs both P and P. To 
simplify the calculations, we will assume that (a) all facets have slopes of ±1, (b) 
facets vanish via binary coarsening, (c) neighboring facet lengths are uncorrelated. 
Additionally, in what follows, it will be helpful to consider Figure [21 which shows 
a representative facet F, and its two neighbors F^ and F^ + \ 
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h- i/2 + L (+) 
h + L/2 
h-L/2 
h + L/2 -l/") 

Figure 2. A diagram illustrating a representative facet and its 
two neighbors. Here h is an arbitrary reference height that occurs 
at the midpoint of the center facet. 




2.1. Flux Law. We begin by observing that, because all slopes are fixed with 
alternating values ±1, then the rate of length change for any single facet is inde- 
pendent of its own vertical velocity, and is instead completely determined by the 
vertical velocity of its immediate neighbors. This general, geometric property can 
be obtained by inspection of Figure O and may be written 

dh i+ i dhi„ 



(2) ^1 = I(_i) 

V ; dt 2 V ' 



dt dt 



where odd (even) facets have negative (positive) slopes. Assume for now that a 
consideration of the facet dynamics ([I]), the existing distribution P, and Equation 
@ allow the derivation of an effective mean facet behavior 

0) m = (§) d), 

which gives the average rate of length change as a function of facet length. Then we 
can use the divergence theorem to write a simple continuity equation for P, letting 
x take the place of L: 

(4) ^ + ^ [P{x)ct>{x)] = 0, 

where 4>(x) now represents the effective lateral velocitey of particles at position x, 
and P(j) then gives the total mass flux. 

2.2. Accuracy: describing coarsening. Since equation (TJ} aims to describe a 
coarsening problem, it must accurately address the primary feature of coarsening 
- the shrinking and vanishing of facets over time. To see if it does, we consider 
Figure [5] again, and imagine that the facet F goes to length 0. We see that three 
things occur: first, the facet F itself vanishes; second, the neighbors of F also 
vanish; and third, a new facet is created which is the merging of the neighbors of 
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F. If these processes are not captured by Eqn. (j4|), then we must add terms to it 
so that it does. 

The first process of facet vanishing is indeed captured by the flux equation Q. 
In the mass-transport framework, vanishing facets are represented by mass particles 
of P(x) flowing out of the domain through the boundary at x = 0. Equation Q 
naturally exhibits this behavior, and allows the easy extraction of the relative rate 
R at which coarsening occurs. This is simply the flux at the origin, given by 

(5) R = -P(O)0(O). 

The second process, neighbor loss, is not captured by To include it, we model 
it as a mass sink 5 1 , which must subtract mass in a way that is probabilistically ac- 
curate. Recalling the assumption that adjacent facet lengths are not correlated, we 
can assume that and -Lj+i are each described by the (probability )distribution 
P(x). Thus, the appropriate sink is 

(6) S(x) = -2P(x). 

The final process of facet creation is also not captured by ((4|. Wc model it, 
in contrast to neighbor loss, as a mass source , which produces facets of length 
Lj_i + Li+i. Since these variables are independent, and each described by P, then 
the sum Li—i + Li + i is described by the joint probability function 

(7) =P 2 (x) = [ P(s)P(x - s)ds, 

Jo 

obtained by integrating a two-point (probability) distribution P(x)P(y) along lines 
of constant x + y. 

The net modifications required by coarsening may now be summed into a single 
term 

(8) C{x) = S{x) + *(ar) = -2P(x) + f P(s)P(x - s) ds. 

Jo 

Since C describes the additional effects of particles leaving the domain (i.e. coars- 
ening), it must be multiplied by R, the rate at which this occurs, and then added 
to Eqn. |(3]). This gives the result 

(9) f = -^p ) + nc(x). 

2.3. Convenience: preserve mass and center of mass. Equation ©, with C 
defined as in ©, now correctly describes the mass redistribution associated with 
coarsening. In particular, it models the two primary features of coarsening - the 
decrease in the number N of facets over time, and the corresponding increase in 
the average length L of those that remain. In the language of a mass-transport 
equation, these correspond repectively to a decrease in the total mass, or zeroeth 
moment M.^°\ and an increase in the center of mass, or first moment Ai^, given 
respectively by 

(o) _ 



(10a) M (V> = / P(x)dx 



(10b) = J^J o xP{x)dx. 
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Indeed, by performing the above integrations against the whole of Eqn. j9]) (rather 
than just against P), we can estabilish that 

(lla) — —J* 

dM (1) f°° 

dm — = l ^ 

These rates of change in A4^ ' and M.^- 1 ' are an accurate description of the relative 
change in the number of facets and average length scale of any initial faceted surface. 
However, the aim of this work is to obtain a description of the steady, dynamically 
scaling distribution of relative facet lengths. This function has not only fixed unit 
mass as P, but also a fixed unit moment. Thus, in this step we perform a dynamic 
non-dimcnsionalization to ensure that both quantities are naturally kept constant. 

To keep X (0) constant, or functionally to make P identically P, we observe in 
Figure [2] that, for each coarsening event, two facets are lost in total. Thus, if the 
rate of coarsening events is R, the net relative rate of mass loss is 2R - indeed, this 
agrees with Eqn. (fTTk ). To keep total mass neutral, we must add mass back into the 
system at the same rate, in a way that is proabilistically accurate. What we want 
to do conceptually is to simply increase the weight of each particle that remains 
in the system, thus preserving total mass without moving any particles. This is 
accomplished mathematically by adding mass particles in precisely the proportion 
that they already exist, as represented by a re- weighting term 

(12) W{x) = 2P(x). 

Finally, to keep A4^ constant, we non-dimensionalize the problem with the 
time- dependent lengthscale L(t) as follows: 

(13) S=^_, 0=^, P = L( t )P, C = L(t)C. 

Since the average lengthscale in the system is constantly increasing due to coars- 
ening, the scaling for x compresses P laterally to keep the center of mass at unity. 
However, to make sure that mass (area) is not also reduced, we also perform the 
scaling shown on P which stretches P vertically, preserving total area, and thus, 
mass. The unknowns <f> and C arc formally scaled as well; <f> is a velocity, and 
so receives the same scaling as x, while C is the sum of two distributions, and so 
receives the same scaling as P. 

Making this substitution and properly applying the chain rulc0 results in a final, 
scale-invariant description of the evolution of the relative length distribution (note 



In particular, application of these scalings to the term dP/dt yields 

dP 1 
~~dt ~* I 



8P 1 dL / dP 
~dt ~ L~dt \ X ~dx~ 



which contains the term (L~ 1 )dL/dt. This is undesireable, as we don't want to keep track of L 
explicitly. However, L <~ A r_1 , and so (L~ 1 )dL/dt = —(N~ 1 )dN/dt, which is a relative rate of 
facet loss. However, we have already defined a similar quantity. Since R describes the relative 
rate of particle loss through the origin, and since two facets vanish for each lost particle (consider 
Figure [2j, 

1 dL _ 1 dN _ 2R 
L dt ~ N dt ~ 
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that S and W cancel) 



(u) ^ = -JL [Pix)(t){x)]+ n 



P 2 (x) + 2(P + x- 
ax 



note that while P was not assumed to necessarily preserve mass nor moment, it 
has been manipulated to do so, and is thus a true probability function. This 
result is quite general, and describes the evolution of the lengthscale distribution 
for any facet dynamics which admits an effective particle velocity (j>. Generally, 
we see that there is a competition between the P2 term, which redistributes mass 
toward infinity, and the (d/dx)(xP) term, which pushes mass toward zero. This 
competition is then modulated by the appropriate dynamics <fi. 

3. Application: Our chosen facet dynamics 

The framework just derived was completely general, describing the coarsening of 
any faceted surface, and indeed any binary system exhibiting binary coarsening. We 
here apply that general framework to the specific facet dynamics ([I]), by deriving 
the appropriate effective flux function 4>(x). Using the dynamics itself, the general 
kinematic form of <j) given in Equation @ , and recalling the diagram in Figure [H 
we perform the following calculation: 

. , dL 1 ( dbS^ dh^ 

(15a) 



dt 2 V dt dt 



(15b) =!((/>(-))_(/>(+) 
(15c) =lfl[2h + L- L^} -±[2h-L + 

(15d) l( L -h L (-) +L M] 



2 V 2 l 

(15e) =IU-^J xP 2 {x)dx 

(15f) 4>(x)=±(x-L(t)) 
(15g) 4>{x) = ^{x-l). 

In step (fTSb). since both and L^ +s> obey the distribution P, the sum L^"' 

and is again modelled by the joint probability function P 2 . The effective 

contribution to <fi is then obtained by performing a weighted integral of possible 
sums x multiplied by their relative prevelance P 2 . Since the total mass of P 2 equals 
unity, that integration describes the center of mass of P 2 , which by considering 
the form of Eqn. (O can be shown to equal 2A4 1 - 1 ' 1 = 2L(t). This result, in turn, 
informs the final nondimcnsionalized result in step {Eg). 

Having calculated 0, we now write the complete evolution equation for the dis- 
tribution of relative facet lengths under the specific facet dynamics ([1]) : 

dP Id 

Remark. We see that the flux p3k ) is negative for facets smaller than the 
average, and positive for facets larger than the average. Broadly speaking, this is 



r x 8 

/ P(s)P(x~s)ds + 2 — (xP(x)) 
Jo dx 
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how coarsening works - small facets shrink and vanish, while large facets grow at 
their expense. Thus, our approximation at least has the right form. 

4. Solution and Comparison with Numerical-Experimental Data 

Solution of Equation (|16p is currently performed numerically, by letting arbi- 
trary initial conditions relax to a steady state; our numerical method is given in 
Appendix [A] This state is unsurprisingly independent of the initial condition cho- 
sen, but surprisingly simple in form - it is simply the exponential distribution 
P(x) = exp(— x), which is easily shown (after the fact) to satisfy Equation (Ti"6| . 
We now proceed to compare the characteristics of this predicted steady state with 
those of the actual steady state found by direct simulation of the dynamics ([1]) ; our 
main results are shown in Figure El 

We begin in Figure by comparing the distribution P itself. The predicted 
exponential distribution is shown in blue; comparison with the green actual distri- 
bution reveals qualitative but not quantitative agreement. In particular, while the 
tail of the predicted distribution is (obviously) exponential, the tail of the actual 
distribution is gaussian. As a consequence, our mean-field steady state exhibits far 
too great an incidence of extremely long facets. 

Seeking the cause of this discrepancy, we next test the accuracy of our effec- 
tive flux function <fi(x) . Figure [3}d shows the contour plot of the distribution of 
length/ velocity pairs p(x,4>). Finding the mean velocity for each length gives the 
statistical </> (dashes), which turns out to compare favorably with the predicted <j) 
(solid). Both are linear with form <f> = a(x — 1), and while the actual slope of 
0.39 differs from the predicted slope of 0.5, they can be made to agree by simply 
scaling time, since every term in Equation <| 1 6j> contains either <p or </>(0). So this 
approximation seems to be valid. 

Finally, we examine the coarsening terms: the sink S(x) in Figure [3fc and the 
source ^(x) in Figure [3ji. These are functionals of P, and so we are not surprised 
that the predicted values (blue) are different from the values calculated from the 
actual steady state (green). However, for both S and ty, we assumed that neighbor- 
ing facet lengths were uncorrclated. If we instead calculate statistically the S and 
^ generated by vanishing facets (that is, the neighbors of vanishing facets), we get 
the curves in red, which are different not only from the predicted quantities, but 
also from the quantites we would have gotten from using the actual distribution P 
and assuming no correlations. 

This suggests that the ultimate culprit is the assumption that neighboring facet 
lengths are uncorrelated. Going back to the simulation, we now measure the corre- 
lation of the lengths of nearby facets as a function of neighbor distance, in Figured] 
There we see a small but significant correlation for at least the first two neighbors. 
This produces the discrepancy between the green and red curves in Figures [3J:,d 
above, and may be responsible for the discrepancy in the tail as well. This result is 
not surprising, as the main weakness of the original LSW theory which inspires this 
work was also a failure to adress correlations; later generalizations which corrected 
this deficiency agreed well with experimental data [40] . 

5. Conclusions 

We have presented a mean-field theory for the evolution of length distribu- 
tions associated with coarsening faceted surfaces. In the spirit of LSW theory, 
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Figure 3. (a) Comparison between the theoretically-predicted 
(blue) and statistically-gathered (green) steady states. The former 
exhibits exponential decay, while the latter is gaussian. (b) Con- 
tour of the statistical distribution of length/ velocity pairs p(x, <j)). 
For each length, mean statistical velocity is plotted as dotted line, 
while the predicted velocity <f>(x) is a solid line. (c,d) Comparison 
of log(— S/2), log(^E') as obtained by various means: from the pre- 
dicted steady state (blue), from the actual steady state (green), 
and from measuring the neighbors of vanishing facets in the actual 
steady state (red). 



a facet- velocity law governing surface evolution is used to establish a characteristic 
length-change law; this in turn leads to a simple continuity equation governing the 
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Figure 4. Statistically-sampled correlation of facet lengths as a 
function of neighbor distance. 



evolution of the probability distribution. However, because the vanishing of any 
facet forces the joining of its two neighbors, this equation must be modified by the 
addition of appropriate terms describing coarsening, including a convolution term 
recalling models of coagulation. Our model therefore serves, apart from the direct 
application to facet dynamics, as a study in the union of these two mechanisms of 
steady statistical behavior. 

The scale-invariant distribution is tracked by introducing an appropriate time- 
dependent non-dimcnsionalization, and the resulting equation is solved by the ex- 
ponential distribution; numerical simulation reveals that any initial condition con- 
verges to this solution. This result unfortunately docs not agree quantitatively 
with the more gaussian distribution obtained by sampling a large surface simulated 
directly under the facet dynamics. Further investigation reveals that the likely cul- 
prit is the assumption that neighboring facet lengths are uncorrelated. Indeed, a 
similar assumption plauged the original LSW theory, and relaxing that assumption 
resulted in much better agreement with experiment. 

However, even with this deficiency, the model captures the essential feature of 
the dynamically scaling state - a particular mass flux law which drives coarsen- 
ing by pushing mass away from unity, moderated by competing terms describing 
coarsening and continuous rescaling, which push mass toward infinity and zero, 
respectively While later improvements to our model addressing neighbor correla- 
tion will undoubtedly increase its predictive capabilities, these same forces will still 
balance in the steady state. The model as presented thus serves as a qualitative 
explanation of the essential features of the scaling state, as well as a guide to further 
reasearch efforts. 
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Appendix A. Numerical Simulation 



This section describes the simulation of Equation (|16| . which is repeated here 
for convenience: 

dP 1 d 



8t 2 X ^- 1 ^ +R 



o 



d 

P(s)P(x -s)ds + 2 — {xP{x)) 
ox 



Since the spatial derivatives here are purely advective, we use an upwinding fi- 
nite difference scheme to approximate them, which easily captures the evolution 
of discontinuous initial conditions, and is accurate to second order in space. Stan- 
dard Simpson's quadrature methods are used for the convolution integral. For 
time-stepping, a Crank-Nicholson integration on the linear terms combined with a 
second-order Adams-Bashforth integration of the integral is quite stable, allwing 
the use of large timesteps to quickly reach the relaxed steady state. However, two 
numerical issues arise that require special mention. 

First, we note that the correct boundary condition on Eqn. (|16p is that P van- 
ishes at infinity. To be able to specify this condition, and also to include the whole 
domain using a finite number of grid points, we perform the change of variables 
£ = cxp(— x), which maps the domain [0, oo) to [1, 0]. The resulting equation, with 
like terms collected, is: 

dP dP f 1 1 

(17) — = [{271 - 0.5) ln(x) - 0.5}x— + {271 - 0.5) P + 71 -P{s)P{£/ s)ds. 
at a£ J ^ s 

Now the entire domain of P is included, and the boundary condition at oo can 
be easily implemented; additionally, gridpoints are distributed unevenly, with a 
bias toward zero to increase resolution of the tail of P. Note, however, that the 
convolution describing P^ now contains an integrablc singularity requiring custom 
quadrature, and additionally requires an interpolation of P to obtain its values at 
the points £/s. 

Second, while Eqn. (|14|) was reached by carefully including terms to conserve the 
mass Ai^ and center of mass A4^\ it turns out that neither is a numerically stable 
quantity. In fact, by performing appropriate multiplications and integrations, we 
can write down a system of ODEs governing the evolution of M.^ and M^: 

(18a) ^ = P(0) (A^> 2 - 1 

(18b) ^=#/M<°'-l 

From this dynamical systems perspective, the point (M.(°\ MS 1 ^) = (1, 1) is indeed 
a fixed point, but an unstable one. This poses no problems for the theory per se, 
as equation (|14[) contains no mechanism for perturbing the solution away from this 
fixed point. However, inevitable numerical error will so perturb any calculation, 
after which such perturbation grows. To avoid this pitfall, it is thus necessary to 
periodically rc-scalc the solution manually back to this fixed point by interpolation. 
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